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Abstract. We present a program called potfit which generates an effective 
atomic interaction potential by matching it to a set of reference data computed 
in first-principles calculations, ft thus allows to perform large-scale atomistic 
simulations of materials with physically justified potentials. We describe the 
fundamental principles behind the program, emphasizing its flexibility in adapting 
to different systems and potential models, while also discussing its limitations. 
The program has been used successfully in creating effective potentials for a 
number of complex intermctallic alloys, notably quasicrystals. 
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1. Introduction 

Classical effective potentials reduce the quantum-mechanical interactions of electrons 
and nuclei in a solid to an effective interaction between atom cores. This greatly 
reduces the computational effort in molecular dynamics (MD) simulations. Whereas 
first principles simulations are limited to a few hundred atoms at most, classical 
MD calculations with many millions of atoms arc routinely performed. Such system 
sizes are possible, because molecular dynamics with short-range interactions scales 
linearly with the number of atoms. Moreover, it can easily be parallelized using a 
geometrical domain decomposition scheme [f , 2], thereby achieving linear scaling also 
in the number of CPUs. 

The study of many problems in materials science and nanotcchnology indeed 
requires simulations of systems with millions of atoms. Quite generally, this is the case 
whenever long-range mechanical stresses are involved. Examples of such problems are 
the study of fracture propagation [3], nano- indentation, or the motion and pinning 
of dislocations. Other problems may be simulated with more moderate numbers of 
atoms, but require very long simulated times, of the order of nanoseconds, an example 
of which is the study of atomic diffusion [4]. In cither case, if large systems and/or 
long time scales are required, classical effective potentials are the only way to make 
molecular dynamics simulations possible. 

The reliability and predictive power of classical MD simulations depend cruicially 
on the quality of the effective potentials employed. In the case of elementary solids, 
such potentials are usually obtained by adjusting a few potential parameters to 
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optimally reproduce a set of reference data, which typically includes a number of 
experimental values like lattice constants, cohesive energies, or elastic constants, 
sometimes supplemented with ab-initio cohesive energies and stresses [5, 6]. In the 
case of more complex systems with a large variety of local environments and many 
potential parameters to be determined, such an approach cannot help, however; there 
is simply not enough reference data available. 

The force matching method [7] provides a way to construct physically justified 
potentials even under such circumstances. The idea is to compute forces and energies 
from first principles for a suitable selection of small reference systems and to adjust 
the parameters of the potential to optimally reproduce them. 

For that purpose, we developed a program called potfit^. By separating the 
process of optimization from the form of the potential, potfit allows for maximal 
flexibility in the choice of potential model and parametrization. 

The underlying algorithms are described in section 2. Section 3 focuses on 
the implementation of the algorithms, followed by details on employing potfit in 
section 4. We discuss advantages and limitations of the force matching method and 
our implementation in section 5, and present our conclusions in the final section 6. 

2. Algorithms 

As mentioned above, potfit consists of two separate parts. The first one implements 
a particular parametrized potential model and calculates from a set of potential 
parameters & the target function that quantifies the deviations of the forces, stresses 
and energies from the reference values. Wrapped around is a second, potential 
independent part which implements a least squares minimization module. As this 
part is completely independent of the potential model and just deals with the list 
of parameters £j, it is fairly straightforward to change the parametrization of the 
potential (tabulated or analytic), or even to switch to a different potential model. 

2. 1 . Optimization 

From a mathematical point of view, force matching is a basic optimization problem: 
There is a set of parameters £j, a set of values bk{£,i) depending on them, and a 
set of reference values &o,fc which the bk have to match. This leads to the well- 
known method of least squares, where one tries to minimize the sum of squares of 
the deviations between the b^ and the 6q fe- I n our case, the reference values can either 
be the components of the force vector /oj acting on each individual atom j, or global 
data Ag^k like stresses, energies, or certain external constraints. We found it helpful 
to measure the relative rather than the absolute deviations from the reference data, 
except for very small reference values. The least squares target function thus becomes 



Z — Zy + Zq 



(1) 



with 





(2) 



j — 1 a—x,y,z 



and 




(3) 
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where Zp represents the contributions of the forces, and Zc that of the global data. 
The (small and positive) Eg impose a lower bound on the denominators, thereby 
avoiding a too accurate fitting of small quantities which are actually not known to 
such a precision. The We are the weights of the different terms. It proves useful for 
the fitting to give the total stresses and the cohesion energies an increased weight, 
although in principle they should be reproduced correctly already from the forces. 
Even if all forces are matched with a small deviation only, those deviations can add 
up in an unfortunate way when determining stresses, thus leading to potentials giving 
wrong elastic constants. Including global quantities in the fit with a sufficiently high 
weight supresses such undesired behaviour of the fitting process. 

As the evaluation of the highly nonlinear target function (1) is computationally 
rather expensive, a careful choice of the minimization method has to be made. We 
chose a combination of a conjugate-gradient-like deterministic algorithm [8] and a 
stochastic simulated annealing algorithm [9]. 

For the deterministic algorithm we take the one described by Powell [8], which 
takes advantage of the form of the target function (which is a sum of squares). By 
re-using data obtained in previous function calls it arrives at the minimum faster 
than standard least squares algorithms. It also does not require any knowledge of the 
gradient of the target function. The algorithm first determines the gradient matrix 
at the starting point in the high-dimensional parameter space by finite differences. 
The gradient matrix is assumed to be slowly varying around the starting point. A 
new optimal search direction towards the minimum is determined by the method of 
conjugate gradients. Then, the target function is minimized along this direction. This 
operation is called line minimization. When the minimum is found, the direction unit 
vector replaces one of the basis vectors spanning the parameter space. The gradient 
matrix is updated only with respect to this new direction, using the finite differences 
calculated in the line minimization. In this way, no finite differences have to be 
calculated explicitly except in the very first step. The line minimization is performed 
by Brent's algorithm [10] in an implementation taken from the GNU Scientific Library 
[11]. 

The algorithm is restarted (including a calculation of the full gradient matrix) 
when either a step has been too large to maintain the assumption of a constant 
gradient matrix, the basis vectors spanning the parameter space become almost 
linearly dependent, or the linear equation involved in Powell's algorithm cannot be 
solved with satisfactory numerical precision. 

The other minimization method implemented is a simulated annealing [12] 
algorithm proposed by Corana [9] . While the deterministic algorithm mentioned above 
will always find the closest local minimum, simulated annealing samples a larger part 
of the parameter space and thus has a chance to end up in a better minimum. The 
price to pay is a computational burdon which can be several orders of magnitude 
larger. 

For the basic Monte Carlo move, we chose adding Gaussian-shaped bumps to the 
potential functions. The bump heights are normally distributed around zero, with 
a standard deviation adjusted so that on average half of the Monte Carlo steps are 
accepted. This assures optimal progress: Neither are too many calculations wasted 
because the changes are too large to be accepted, nor are the steps too small to make 
rapid progress. 
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2. 2. Potential models and parametrizations 

The simplest effective potential is a pair potential, which only depends on interatomic 
distances. It takes the form 

N 

V^=E<W^), (4) 

i,j<i 

where is the distance between atoms i and j, and (f> SiSj is a potential function 
depending on the two atom types Sj and Sj. This function can cither be given in 
analytic form, using a small number of free parameters, like for a Lennard- Jones 
potential, or in tabulated form together with an interpolation scheme for distances 
between the tabulation points. Whereas the parameters of an analytic potential can 
often be given a physical meaning, such an interpretation is usually not possible 
for tabulated potentials. On the other hand, an inappropriate form of an analytic 
potential may severely constrain the optimization, leading to a poor fit. For this 
reason, we chose the functions (f> to be defined by tabulated values and spline 
interpolation, thus avoiding any bias introduced by an analytic potential. This choice 
results in a relatively high number of potential parameters, compared to an analytic 
description of the potentials. This is not too big a problem, however. Force matching 
provides enough reference data to fit even a large number of parameters. The potential 
functions <j) only need to be defined at pair distances r between a minimal distance 
r m ; n and a cutoff radius r cut , where the function should go to zero smoothly. 

We found pair potentials to be insufficient for the simulation of complex metallic 
alloys. More suited are EAM (Embedded Atom Method [13, 5]) potentials, also 
known as glue potentials [14], which have many advantages over pair potentials in 
the description of metals [14]. 

EAM potentials include a many-body term depending on a local density nf. 

V = ^ (j>8 iSj (rjj ) + J~] U Si (nf) with n, = ^/) S] (r,j). (5) 

i,j<i i j^i 

Hi is a sum of contributions from the neighbours through a transfer function p s , , and 
U Si is the embedding function that yields the energy associated with placing atom i at 
a density n^. Again, all functions are specified by their values at a number of sampling 
points. 

The parameters £j specifying a tabulated potential are naturally the values at the 
sampling points. Due to the nature of spline interpolation, either the gradient or the 
curvature at the exterior sampling points of each function can also be chosen freely. 
Depending on the type of potential one can keep the gradients fixed, or adapt them 
dynamically by adding them to the set of parameters 

The EAM potential described by (5) has two gauge degrees of freedom, i.e., two 
sets of parameter changes which do not alter the physics of the potential: 

p s (r) -> np s (r), 

and 

^s iSj (r) -> <f>s iSj {r) + A Si p Sj (r) + X S] p Si (r), 

U s ,( n i) ~* U Sz (ni) - \ Si Tii. 
According to (6), the units of the density rij can be chosen arbitrarily. We use this 
degree of freedom to set the units such that the densities rij computed for the reference 
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configurations arc contained in the interval (—1; 1], but not in any significantly smaller 
interval. The transformation (7) states that certain energy contributions can be moved 
freely between the pair and the embedding term. An embedding function U which is 
linear in the density n can be gauged away completely This also makes any separate 
interpretation of the pair potential part and the embedding term void; the two must 
only be judged together. The latter degeneracy is usually lifted by choosing the 
gradients of the Ui(ni) to vanish at the average density for each atom type, potfit also 
uses this convention when exporting potentials for plotting and MD simulation. As 
the average density might change during minimization, potfit internally uses a slightly 
different gauge: It requires that the gradient vanishes at the center of the domain of 
the respective embedding function. 

potfit can perform the transformations (6,7) periodically on its own, thus 
eUminiating the need to fix the gauge by an additional term in the target function 
(1). Unfortunately, for tabulated functions the transformations cannot be performed 
exactly due to the nature of spline interpolation. A change of gauge therefore can 
lead to an increase of the target function, which is why we suppress such gauge 
transformations in the very late stages of a minimization. 

3. Implementation 

potfit is implemented in ANSI C. While the user may specify most options in a 
parameter file read when running the program, some fundamental choices must be 
made at compile time, like for example the potential model used, or whether to allow 
for automatic gauge transformations in EAM potentials. This is a compromise between 
convenience and computation speed. Compile time options can be selected by passing 
them to the make command, and thus do not require any changes of the source files. 
For solving the linear equations in Powell's minimization algorithm, potfit makes use of 
routines from the LAPACK library [15], which must be installed separately, probably 
together with the BLAS library [16] LAPACK is based on. 

3.1. Parallelization and optimization 

The program spends almost all CPU time in calculating the forces for a given potential; 
finding a new potential to be tested against the reference data takes only a tiny 
fraction of that time. Thus, the only way to improve performance is to reduce the 
total time needed for the force computations, either by minimizing their number, or by 
making each force computation faster. Powell's algorithm leaves only little room for 
further reduction of the number of force evaluations. One could for instance adjust the 
precision required in a line minimization. If the tolerance is too small, time is wasted 
in refining a minimum beyond need, whereas an insufficent precision may stop too far 
from the minimum, thus requiring more steps in total. The choice of this tolerance 
was made empirically. 

Much more time can be saved by parallelizing the calculation of forces, energies, 
and stresses for a given potential. This is done in a straightforward way: As 
the forces, energies, and stresses of the different reference configurations can be 
computed independently, we simply distribute the reference configurations on several 
processes. Before the force computation, the potential parameters are distributed to all 
processes, and afterwards the computed forces, energies and stresses are collected. The 
communication is performed using the standard Message Passing Interface (MPI [17]). 
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This simple parallclization scheme works well as long as the number of configurations 
per process does not drop below 10 to 15. Otherwise, the communication overhead 
starts to show up, and load balancing problems may appear. A shared memory 
OpenMP parallclization also exists, but produces inferior results. 

In force matching, the reference configurations stay fixed. Therefore, all distances 
between atoms remain fixed, and ■potfit can use neighbour lists, which need to be 
computed only once at startup. In fact, for each neighbour pair all data required 
for spline interpolation are pre-computed, allowing for a fast lookup of the tabulated 
functions. This data needs to be recomputed only when the tabulation points of a 
function are changed. 

3. 2. Input and output files 

Tabulated potential functions can be specified with equidistant or with arbitrary 
tabulation points. For equidistant tabulation points, the boundaries of the domain 
and the number of sampling points of each function are read from the potential file, 
followed by a list of function values at the sampling points and the gradients at the 
domain boundaries. In the case of free tabulation points, only their number is specified 
at the beginning of the potential file, followed by a list of argument-value pairs and 
again the gradients of the potential functions at the domain boundaries. 

Reference configuration files contain the number of atoms, the box vectors, the 
cohesive energy, and the stresses on the unit cell, followed by a list of atoms, with 
atom species, position and reference force for each atom. Such reference configuration 
files can simply be concatenated. 

potfit was designed to cooperate closely with the first-principles code VASP 
[18, 19] and with IMD [20], our own classical MD code. VASP, which is a plane 
wave code implementing ultrasoft pseudopotentials and the Projector- Augmented 
Wave (PAW) method [21, 22], is used to compute the reference data for the force 
matching, whereas the resulting potentials are intended to be used with IMD. For this 
reason, potfit provides import and export filters for potentials and configurations to 
communicate with these programs. These filters are implemented as scripts, which 
can easily be modified to interface with other programs. 

4. Results and validation 

As a first test, potfit should be able to recover a classical potential from reference data 
computed with that potential. For this test, we used snapshots from several molecular 
dynamics runs as reference structures, first for a Lennard- Jones fee solid, then for a 
complex Ni-Al alloy simulated with EAM potentials [23]. In order to ensure that all 
reference data presented to potfit is consistent, the potentials were approximated by 
cubic spline polynomials, in the same way as potfit represents the potentials. With 
such reference data and starting with vanishing potential functions, potfit could in both 
cases perfectly recover the potentials. This test therefore demonstrates the correctness 
of the program. One should keep in mind, however, that reference data from ab-inito 
computations often cannot be reproduced perfectly by any classical potential. 

Our primary research interest are quasicrystals [24] and other complex metal 
alloys, for which good potentials are hardly available, potfit has been developed in 
order to generate effective potentials for such complex metal alloys, which feature large 
(or even infinite) unit cells, several atom species, and a wide variety of different local 
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environments. So far, force matching had been used mainly to determine potentials 
for monoatomic metals and a small selection of relatively simple binary alloys. 

As a first application beyond simple alloys, we have developed potentials for the 
quasicrystalline and nearby crystalline phases in the systems Al-Ni-Co, Ca-Cd, and 
Mg-Zn. Due to the complexity of the structures and also due to the choice of tabulated 
potential functions, a relatively large number of potential parameters is required. This 
is especially true for ternary EAM potentials, which comprise 12 tabulated functions, 
with 10-15 tabulation points each. Correspondingly, a relatively large amount of 
reference data is required. A computationally efficient implementation of the force 
matching method is therefore essential. It turned out that potfit scales well under 
those circumstances and is up to its task. 

Although the potentials to be generated are intended for (aperiodic) quasicrystals 
and crystals with large unit cells, all reference structures have to be periodic crystals 
with unit cell sizes suitable for the ab-initio computation of the reference data. On 
the other hand, the reference structures should approximate the quasicrystal in the 
sense, that all their unit cells together accommodate all relevant structural motifs. 
To do so, they must be large enough. For instance, the quasicrystalline and related 
crystalline phases of Ca-Cd and Mg-Zn consist of packings of large icosahcdral clusters 
in different arrangements. Reference structures must be able to accomodate such 
clusters. A further constraint is, that the unit cell diameter must be larger than the 
range of potentials. We found that reference structures with 80-200 atoms represent 
a good compromise between these requirements. 

Starting from a selection of basic reference structures, further ones were obtained 
by taking snapshots of MD simulations with model potentials at various temperatures 
and pressures. Also samples which were strained in different ways were included. 
For all these reference structures, the ab-initio forces, stresses and energies werde 
determined with VASP, and a potential was fitted to reproduce these data. As 
reference energy, the cohesive energy was used, i.e., the energies of the constituent 
atoms was subtracted from the VASP energies. Instead of absolute cohesive energies 
one can also use the energy relative to some reference structure. Once a first version of 
the fitted potential was available, the MD snapshots were replaced or complemented 
with better ones obtained with the new potential, and the procedure was iterated. 

As expected, no potential could be found which would reproduce the reference 
data exactly. During the optimization, the target function (1) does not converge 
to zero, which indicates that quantum mechanical reality (taking density functional 
theory as reality) is not represented perfectly by the potential model used. The forces 
computed from the optimal potential typically differ by about 10% from the reference 
forces, which seems acceptable. For the energies and stresses a much higher agreement 
could be reached. Cohesion energy differences for instance can be reproduced with an 
accuracy better than 1%. 

The generated potentials were then used in molecular dynamics simulations to 
determine various material properties, such as the melting temperature and the elastic 
constants, for which values consistent with experiment were obtained. The Ca-Cd 
potentials were especially tuned towards ground-state like structures, whose energies 
are reproduced with high accuracy, in agreement with ab-initio results. Details of 
these applications can be found in [4, 25, 26, 27]. Probably the best tested EAM 
potential constructed with potfit was obtained by Rosch, Trcbin and Gumbsch [28]. 
This potential is intended for the simulation of crack propagation in the C15 Laves 
Phase of NbCr2, and has undergone a broad validation. These authors calculated the 
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lattice constant, the clastic constants and the melting temperature and compared 
these values to experimental and ab-initio results with reasonable success. They 
also studied relaxation of surface atoms, surface energy and the crack propagation 
in NbCr2. According to the authors [28], the force- matched potentials created with 
potfit clearly outperform previously published potentials. But this example also shows 
[28] that a large number of fitting-validation cycles arc usually required, before a usable 
and satisfactory potential is obtained. This makes force matching a time-consuming 
and tedious process. 

5. Discussion 

5. 1 . Transferability 

It should be kept in mind that force-matched potentials will only work well in 
situations they have been trained to. Therefore, all local environments that might 
occur in the simulation should also be present in the set of reference configurations. 
Otherwise the results may not be reliable. Using a very broad selection of reference 
configurations will make the potential more transferable, making it usable for many 
different situations, e.g. for different phases of a given alloy. On the other hand, 
giving up some transferability may lead to a higher precision in special situations. By 
carefully constraining the variety of reference structures one may generate a potential 
that is much more precise in a specific situation than a general purpose potential, 
which was trained on a broader set of reference structures. The latter potential, on the 
other hand, will be more versatile, but less accurate on average. Finding sufficiently 
many suitable reference structures might not always be trivial. For certain complex 
structure like quasicrystals, there may be only very few (if any) approximating periodic 
structures with small enough unit cells. 

5. 2. Optimal number and location of sampling points 

Each reference database has an optimum number of parameters it can support. Using 
too few parameters, the potential functions lack flexibility. On the other hand, 
exceeding this number may lead to overfitting beyond the limit of the potential 
model, potfit cannot determine that optimal number automatically, but there is a 
simple strategy the user can employ. The set of reference configurations is split in 
two subsets, one of which is used for fitting and the other for testing the potential. 
If the root-mean-square (rms) deviation of the test set significantly exceeds that of 
the fitting set, the database is probably ovcrfittcd [29]. By starting with a relatively 
low number of parameters, that is increased as long as the rms of the testing stage 
decreases, one can arrive at the optimal number of parameters [30]. 

This strategy also helps in dealing with oscillatory artefacts of the spline 
interpolation: If the sampling points are not spaced too densely, and there is enough 
data to support each tabulation point, artificial wiggles are suppressed, potfit provides 
the frequency with which each tabulation interval is accessed during an evaluation of 
the target function (1). With this information, sampling point density can be reduced 
for distances that do not appear frequently enough in the reference configurations. 
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5. 3. Number of atom types and choice of reference structures 

The most obvious impact of an increasing number of atom types is the corresponding 
increase in the number of potential parameters. For instance, an EAM potential 
for n atom types requires n(n + l)/2 + 2n tabulated functions, each with 10 to 15 
tabulation points. For a ternary system, this already amounts to the order of 150 
potential parameters. Whereas such a number of parameters can still be handled, an 
increasing number of atom types leads to yet another problem, which is more serious. 
To see this, it must be kept in mind that any potential function depending on the 
interatomic distance must be determined for the entire argument range between r m j n 
and r cut . If tabulated functions are used, for each tabulation interval there must be 
distances actually occuring in the reference structures, for otherwise there are potential 
parameters which do not affect the target function, and which conseqently cannot be 
determined in the fit. The requirement that all distances for all combinations of atom 
types actually occur in the reference structures becomes especially problematic if the 
atoms of one type form only a small minority, in which case some distances between 
such atoms might be completely absent in all reasonable reference structures. If the 
number of atom types is large, there is unfortunately always at least one element 
which is a minority constituent. In such situations it might be unavoidable to use a 
much broader selection of reference structures with varying stoichiomctry, instead of a 
fixed stoichiomctry with a minority constituent. It might even be necessary to include 
energetically less favourable configurations to provide a complete set of reference data. 

Another solution would be to use a non-local (or less local) paramctrisation of 
the potential functions, like a superposition of broad gaussians or functions given by 
analytic formulae. Changing one parameter can then affect the function over a broader 
range of arguments, making it again possible to fit the function even if only sparse 
information on it is provided by the reference data. Potentials represented in this way 
would also not suffer from the wiggle artefacts of spline interpolation described above. 

5.4- Experimental values as reference data 

potfit does currently not use experimental data during force matching. The potentials 
are determined exclusively from ab-initio data, which means they cannot exceed the 
accuracy of the first principles calculations. While it is possible, in principle, to support 
also the comparison to experimental values, we decided against such an addition. 
For once, available experimental values can often not be calculated directly from the 
potentials, so determining them would considerably slow down the target function (1) 
evaluation. Secondly, experimental values often also depend on the exact structure 
of the system, which in most cases is not completely known beforehand for complex 
structures, for instance due to fractional occupancies in the experimentally determined 
structure model. A better way to use experimental data is to test whether the newly 
generated potentials lead to structures that under MD simulation show the behaviour 
known from experiment. 

6. Conclusion 

Large scale molecular dynamics simulations are possible only with classical effective 
potentials, but for many complex systems physically justified potentials do not exist so 
far. Our program potfit allows the generation of effective potentials even for complex 
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binary and ternary intcrmctallics, adjusting them to ab-initio determined reference 
data using the force matching method. Potentials for several complex intermetallic 
compounds have been generated, and were successfully used in molecular dynamics 
studies of various properties [4, 25, 26, 28]. It should be emphasized, however, that 
constructing potentials is still tedious and time-consuming. Potentials have to be 
thoroughly tested against quantities not included in the fit. In this process, candiate 
potentials often need to be rejected or refined. Many iterations of the fitting-validation 
cycle are usually required. It takes experience and skill to decide when a potential is 
finished and ready to be used for production, and for which conditions and systems 
it is suitable, potfit is only a tool that assists in this process. Flexibility and easy 
extensibility was one of the main design goals of potfit. While at present only pair 
and EAM potentials with tabulated potential functions are implemented in potfit, it 
would be easy to complement these by other potential models, or to add support for 
differently represented potential functions. 
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